Elastic scattering theory and transport in graphene 
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Electron properties of graphene are described in terms of Dirac fermions. Here we thoroughly 
outline the elastic scattering theory for the two-dimensional massive Dirac fermions in the presence 
of an axially symmetric potential. While the massless limit is relevant for pristine graphene, keeping 
finite mass allows for generalizations onto situations with broken symmetry between the two sublat- 
tices, and provides a link to the scattering theory of electrons in a parabolic band. We demonstrate 
that the Dirac theory requires short-distance regularization for potentials which are more singular 
than 1/r. The formalism is then applied to scattering off a smooth short-ranged potential. Next we 
consider the Coulomb potential scattering, where the Dirac theory is consistent for a point scatterer 
only for the effective impurity strength below 1/2. From the scattering phase shifts we obtain the 
exact Coulomb transport cross-section in terms of the impurity strength. The results are relevant 
for transport in graphene in the presence of impurities that do not induce scattering between the 
Dirac points in the Brillouin zone. 

PACS numbers: 81.05.Uw 72.10.-d 73.63.-b 73.40.-c 



INTRODUCTION AND OUTLINE 

Graphene, a layer of Carbon atoms arranged in a hon- 
eycomb lattice, has been long known for its peculiar 
electronic dispersion, equivalent to that of massless two- 
dimensional (2D) Dirac fermions. This system was first 
considered in tight-binding approximation by Wallace 1 
and by McClure. 2 For a long time, graphene mono- 
layer served as a low-dimensional toy model where Dirac 
fermions appear naturally. 4,5 Significant interest to this 
material arose in 1990s fueled by the discovery of carbon 
nanotubes. 3 The field has experienced an even stronger 
surge of interest since 2004, when graphene monolayers 
were obtained experimentally. 6 The outstanding quality 
of graphene monolayers and few-layered samples is mani- 
fest in high mobility resulting in ballistic conductance on 
fxm scale, and in quantized Hall effect. 6-10 

Recent electron transport measurements 7-9 show that 
the mobility in graphene is approximately independent 
of the carrier density (i.e. conductivity grows propor- 
tional to the density). The effects of various kinds of 
the potential disorder on transport in graphene have 
been considered in a number of works. 11-27 Arguably, 
the density-independent mobility originates mainly due 
to the Coulomb impurities in the substrate. 15-17 

The importance of the smooth potential disorder for 
the transport in graphene prompts the development of 
the scattering theory for the 2D Dirac fermions, both 
massless and massive. Physically, a nonzero mass can 
arise due to an external perturbation that distinguishes 
between the sublattices; recent ab-initio density func- 
tional calculations predict the Dirac gap of 53 meV when 
placing graphene monolayer on a hexagonal boron ni- 
tride substrate. 28 Somewhat similar perturbation occurs 
in bilayer graphene, 29-31 although in this example the 
spectrum is not exactly of the Dirac form. An interest- 



ing possibility for the Dirac gap opening up due to the 
spin-orbit coupling was considered by Kane and Mele. 32 

The purpose of this work is to thoroughly outline the 
elastic scattering theory for 2D Dirac fermions in the 
axially-symmetric potentials. Such a 2D formalism is 
built essentially following Ref. 33, the classic reference 
for scattering in 3D Dirac systems. The connection to 
the transport in graphene in the presence of potential 
centers whose field is smooth on the lattice scale is es- 
tablished via the transport cross section. 

We start from the basic facts about the 2D Dirac 
model, and write the normalized spinor plane and spher- 
ical waves (Sec. I). In Sec. II we study the properties of 
the radial solutions for the 2D Dirac spinors, define the 
scattering phase shifts, and link them to the differential 
and transport cross sections. We also derive the Born 
approximation for the 2D Dirac spinors, as well as out- 
line analytical properties of the radial solutions on the 
complex energy plane. 

One of the observations made in Sec. II is that the 
Dirac problem, both massive and massless, requires a 
short-distance regularization whenever the external po- 
tential is more singular than 1/r. Classically, this cor- 
responds to falling into the potential center. For such 
singular potentials the purely Dirac formalism is inappli- 
cable, and the lattice scale physics starts playing a role. 

In Sec. Ill we consider scattering off a potential local- 
ized within a finite radius that exceeds the lattice con- 
stant but is smaller than the particle wavelength. 

In Sec. IV we focus on one of the most important scat- 
tering problems for graphene — that for the potential 
U = —hvx a jr. The latter problem has so far been 
treated in the Born approximation. 15-17 The exact solu- 
tion is presented in detail for both massless and massive 
cases, and for both signs of the impurity charge a. The 
asymptotic behavior of the scattering solutions and scat- 
tering phase shifts are studied, with the particular atten- 
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tion paid to the "ultrarelativistic" and "nonrelativistic" 
limits, relevant correspondingly for the pristine graphene, 
and for the graphene layer with broken sublattice sym- 
metry, or for the electrons in a semiconductor with a 
parabolic band. 

The 1/r problem deserves a special consideration, as 
it is a borderline case for the falling into the potential 
center. It has been known 34 that the solutions in the 1/r 
potential become singular whenever the "fine structure 
constant" a > a c — 1/2, i.e. at even smaller value than 
that in 3D (a c = 1). In this respect the physics of the 
Coulomb impurity in graphene, depending on the dielec- 
tric environment, may correspond to the "supercritical" 
relativistic heavy atom (that with Z > 137). 35-38 

Finally, In Sec. V we use the exact scattering phases 
to calculate the transport cross section for the subcriti- 
cal Coulomb impurity, and compare the exact result with 
the Born approximation. Our main finding there is that, 
for a given carrier charge, the attracting impurity scat- 
ters more effectively than the repelling one. This should 
be contrasted with known exact 2D and 3D nonrelativis- 
tic scattering results in a 1/r potential, where such an 
asymmetry does not take place. 



I. FREE ELECTRONS IN GRAPHENE 
A. Model 

The key features of electron dispersion in an ideal 
graphene monolayer can be summarized as follows. 3 With 
the two cites per unit cell, graphene's 7r-electron band 
has the two inequivalent points in the Brillouin zone, at 
which the electron and hole subbands just barely touch. 
At these so-called Dirac points K and K' , the carrier dis- 
persion is linear and electron-hole symmetric, e(p) oc ±p. 
Separately at the K and K' points, the wave function in 
the effective-mass approximation has a spinor structure, 
its two components corresponding to the two sublattices. 
This spinor obeys the massless Dirac equation. The low- 
energy states near the K and K' points are decoupled in 
a pristine graphene monolayer. Provided the material is 
subject to external fields that are adiabatic on the lat- 
tice scale, the low-energy properties can be understood in 
terms of the Nf = 2 sp i„ x 2 va n ey = 4 independent Dirac 
fermion polarizations. 

Near a given Dirac point, the free effective- mass Hamil- 
tonian 

Ho = -ihv(T 3 aid x + <T 2 d y ) + Acr 3 , (1.1) 

where v « 1 x 10 6 m/s is the graphene Fermi velocity, 
°X2,3 are the Pauli matrices that act in the spinor space 
corresponding to the two sublattices of a honeycomb lat- 
tice, while T3 = ±1 distinguishes between the K and K' 
Dirac points. Everywhere in this work we consider the 
dynamics on the scale much larger than the graphene's 
lattice constant and neglect scattering between the Dirac 



points. Hence, without loss of generality, we set T3 = +1 
in what follows. 

In the Hamiltonian (1.1) we also introduced the gap 
(the Dirac mass) 

A = Mv 2 . (1.2) 

Although this term is absent by symmetry in an ideal 
graphene monolayer, it can become imporant when the 
symmetry is reduced. Without loss of generality here we 
set M > (working at zero magnetic field, we are not 
concerned with the parity anomaly effects ' 5 ). 

The Hamiltonian (1.1) corresponds to the Lagrangian 
[h = v = l,T 3 = +1] 

£o = ^(H^-M)V, $ = 1>h°, (1-3) 

where the Dirac matrices 7 = 0-3, 7 1 — io<i and 
7 2 = -icru such that {7^, 7"} = 2g» v with g^ u = 
diag (1, — 1, — 1). In this "relativistic" notation the 
Lorentz-invariant fermion current is 

= ^ = (P, J) , (1.4) 

where the number density and the current are 

p = VV, 3 = ^crip, (1.5) 

with cr = (<7i, 02)- 

B. Spinor plane waves 

Consider the eigenproblem Hoip — tip, Eq. (1.1), where 
we represent the two-component spinor 

V>(r) - Q . (1.6) 

The components of ip satisfy [we set fi = v = 1 in what 
follows] 

( M P X ~ i Py \ (<P\_ f (<P\ (I -N 

{px+iPv -M J [xj W ' 

where p x = —id x and p y = —id y . In the plane wave 
basis (tp x) T — Ue, P e lpr the differential operators become 
components of the momentum eigenvalue p, yielding the 
relativistic dispersion 

e = ±vV + M2, P=^pl+V 2 y , (1-8) 

where ± distinguishes between the particle and hole sec- 
tors. The conventional normalization "one particle in 
a unit volume" 33 J M = p^/e = (l,v p J, where v p = 
de/dp = p/e is the velocity, requires Ve.p^e.p = M/e 
or, equivalently, V^pV^p — 1j yielding 

4, -u u+u --^( ^FTiWT \ 

^ e;P e;P ' M±H;P ~ y/2\T\ {±VW^M\e^J- 

(1.9) 
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Here 9 p = arg (p x + ip y ), and the upper and lower signs 
refer to the electron (e > M) and hole (e < —M) parts of 
the spectrum, ± = signe. The absolute values under the 
square roots are introduced to describe both sectors. The 
factor w = is an overall phase that has a meaning of 
the "nonrelativistic" particle's wave function in the rest 
frame, e = M. 



C. Spinor spherical waves 

For the purpose of developing the scattering theory, 
below we introduce the spherical wave basis of eigenstates 
of the problem (1.1). 

First recall that the 2D nonrelativistic scalar particle 
with fixed absolute value of the momentum p = |p| and 
fixed projection m = 0, ±1, ±2, ... of angular momentum 
on the z-axis (perpendicular to the plane) is described 
by the spherical wave 



V pm {r) = § m {6)R pm {r) . 
Here the angular harmonics 



2tt 



to = 0,±1,±2, ... , 



(1.10) 



(1.11) 



and the radial functions R pm (r) satisfy the radial 
Schrodingcr equation 

1 d ( d \ m 2 2 
~rdr~ \d^ Rpm ) + — Rpm = P Rpm > (L12) 

which reduces to the Bessel equation 

p 2 R'; p + P R' p + (p 2 -TO 2 )i? = 0, p = pr. (1.13) 

The solutions are the Bessel functions J m {p) and the 
Neumann functions Y m (p), whose asymptotic behavior 



nr 
V n p 
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(>- 


m7r 
~ 


nr 
V ftp 


sin | 







Y m (p > m) 

Their corresponding short-distance behavior is 

Jm(pr) ^i(y) m 



and 



Y m (pr) 



T(m) /_2_V 
7T \pr J 



to > 0; 



Mn( lE pr/2), to = 0, 



(1.15a) 



(1.15b) 



where \nj E ~ 0.577... is the Euler's constant. [For 

to < 0, use J_ m = (-) m J m and Y_ m = {-) m Y m ] 
Thus for the wave regular at r = 0, one chooses R m cx 



J m (pr), normalized according to J^rdr R pm R P >m' — 
2n6 mm ,6(p-p'): 



R v 



r- T . . 2 / TO7T 7T\ 

/27rp x J m {pr) ~ —j= cos [pr - — - -J . 



(1.16) 

Coming back to the relativistic case, one notes that 
both the isospin \<j and the angular momentum l z — 
—id$ = ~i(xd y — yd x ) do not commute with the Hamil- 
tonian (1.1): 

[l z , Ho] = i<r x p , [\a z , H ] = -i(r x p , 

thus a state cannot be characterized by their values. (In 
fact, the spherical spinor components will have different 
eigenvalues of l z .) The conserved quantity is the "isospin- 
orbital" momentum around the z-axis, 11 



J 



(1.17) 



Also, similar to the 3D case, 33 the parity of a state is 
conserved: Under inversion r ^ — r [i.e. 6 ^> 9 + tt for 
the polar angle, and tp — > 7°^], the spinor components 
(1.6) transform as <p(r) — > <p(— r) and x( r ) — * — r )- 
The spinor Vpm w hl have the definite parity (— ) m if its 
components 



Vv»(r) = 



^G pm (r)$, 



(1.18) 



have the angular parts correspondingly with l z — m and 
i 2 = to + 1. The factor of z here is chosen for later 
convenience. 

Consider now the radial parts F pm (r) and G pm (r) as- 
suming |e| > M. From the relations (1.7), it follows that 



1 d 
r dr 
ld_ ( dG 

r dr 



dF„ 



dr 



r 2 P m 



p 2 F pmi (1.19a) 



dr 



(m+ If 



G pm = P 2 G pm . (1.19b) 



The equations (1.19) are of the radial Schrodingcr 
form, Eq. (1.12). Thus F pm = AR pm (r) and G pm = 
BR P! m+i(r). For the spinor wave regular at the ori- 
gin, one chooses the radial functions in the form (1.16). 
To find A and B, we consider the limit pr — > oo, 
when the wave function is approximately a plane wave 
in the direction of f . Using the asymptotic behavior 
(1.16) and the relations (1.9) between the components 
of the plane wave, find B/A = ± yfe - M|/|e + M\, 
where ± = signe. Requiring the overall normalization 
JcPrtpp^p,^ = 2Tt8 mm >5(p - p'), obtain the spinor 
spherical wave 



VVm( r ) 



^\e + M\R pm (r)f> m (0) N 
±iy/\e - M\Rp im+1 (r)$ m+ i {9), 



(1.20) 
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whose parity is (— ) m . Here R 
the spinor regular at r = 
for the spinor singular at r - 



prn, 

R 

o'. 



(r) = y/2n pJ m (pr) for 
pm (r) = yJ2irpY m (pr) 
and same applies for 



Rp,m+i(r)- The spinor (1.20) is also an cigcnstate of the 
operator (1.17) with the eigenvalue 



(1.21) 



Below we will often use the eigenvalue (1.21) instead of 
the orbital number m to label states, eigenvalues or phase 
shifts; e.g. for the spinors (1.20), ip pm (r) = i[) p j(r). 



II. POTENTIAL SCATTERING 

A. Equations for the radial functions 

We consider elastic scattering off an axially-symmetric 
external scalar potential U(r). The Hamiltonian 

H = -i(a 1 d x +a 2 d y )+M<j 3 + U(r). (2.1) 

In the spinor components (1.6), Eq. (2.1) reads 



(e - M - U)ip = (p x - ip y )x 
(e + M - U)x = (p x + ip v )v , 



(2.2) 



where (p x , p y ) = (—id x , —id y ) are differential operators. 

The crucial symmetry of the problem (2.1) is the con- 
servation of the total orbital momentum (1.17), since 



[l x , H] = -[\a z , H}=i<TXp 



[j z , H] = (2.3) 



for any axially-symmetric U(r). This property allows us 
to work in the spherical basis of the form (1.18). Taking 
the spinor (1.6) in the form (1.18), and using <& m+ i = 
e j9 $ m , and p x ± i Py = e ±l6 (-id r ± ±d g ), where 9 = 
arg(x + iy), obtain the following equations for the radial 
functions F and G: 



If - -F + (e + M - U)G 

dr r 



dG m + 1 
— — I G 

dr r 



(e-M- U)F 



0, 
0. 



(2.4a) 
(2.4b) 



Everywhere here it is implied that the functions F and 
G correspond to the angular momentum (1.21), e.g. F = 
F m = Fj; the index j or m will be often suppressed 
for brevity. Eqs. (2.4) have been derived by DiVincenzo 
and Mele 11 for the case M = 0. In the absence of the 
potential, U = 0, Eqs. (2.4) are equivalent to Eqs. (1.19). 

It is often convenient to represent Eqs. (2.4) in a more 
symmetric form, using the eigenvalue (1.21), 

(Fy/r)' r - -{Fy/r) + (e + Af - U)(Gy/r) = 0,(2.5a) 
r 

{Gy/f)' r + ^{(Gy/r)-{e-M-U){Fy/r) = 0.(2.5b) 

Eqs. (2.4) and (2.5) are valid both for the continuous 
and for the discrete spectrum (present for M ^ 0). In 



the massless limit M — > 0, Eqs. (2.4) and (2.5) acquire 
the following symmetry: For any U(r), if a pair (F, G) is 
the solution for a given j, then the pair (G, —F) is the 
corresponding solution for j — ► — j, i.e. 



M = 



(2.6) 



This symmetry is also present for the 3D massless Dirac 
fermions [Ref. 33, Sec. 38]. 

The asymmetry between the spinor components associ- 
ated with the finite Dirac mass M is stressed by rescaling 
the radial functions in accord with Eq. (1.20), 



lei > M : 



Fy/r = y/±(e + M)F, 



Gyfr~ = ±y/±(e-M)G, 
where ± is signe. Then Eqs. (2.5) take the form 



K- 1 f + p 



g; + -g 

r 



P 



U 



e + M 

U 
e-M 



G 
F 



0, 
0. 



(2.7) 

(2.8a) 
(2.8b) 



Here p = y/e 2 — M 2 . The same can be done for the 
discrete spectrum, 



\e\ < M 



Fy/r~= y/M + eF, 



(2.9) 



Gy/r = yfW~~€ G . 

Introducing A = y/M 2 — e 2 , the corresponding equations 



F' - -F + X 

r 

G' + -G + A 

r 



U 



1 + 



M + e 
U 



M - 



G = 0, 
F = 0. 



(2.10a) 
(2.10b) 



Finally, we reduce the system of first-order equations, 
say Eqs. (2.8), to an equivalent second-order equation. 
The latter can be written either for F or for G, as follows: 



F" + 



U' 



F' 



'-F 



e + M -U 



(e - U) 2 - M 2 + 



(2.11) 

[the corresponding equation for G would have j — > —j] . 
Eq. (2.11) is reduced to the familiar Schrodinger form 



*" + 2[£-F(r)]* = 0, E = p 2 /2, 



(2.12) 



via the substitution F = y/e + M — U ^. Similar to the 
3D case, the potential V = V\ + V2 splits into the Klein- 
Gordon part V\ , and the part V2 responsible for the Dirac 
"spin" effects: 



1 i 2 



V, - *U{r)--U'+ 2r2 



J 



(2.13) 



V 2 = 



U" 



e + M -U 2\e + M -U 



U' 



2W 



e + M -U 
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It should be clear that the "spin" -orbit coupling coming 
from the potential Vi has nothing to do cither with the 
real SU(2) spin of electrons, or with having the two Dirac 
points. Rather, it is a consequence of a two-component 
spinor structure of the electron wave function due to the 
existence of the two sublattices in a honeycomb lattice. 
We also note that the Schrodinger "energy variable" E in 
Eq. (2.12) has the dimension of [energy] 2 , same as that 
of the potential V. In effect, Eq. (2.12) is a square of the 
original Dirac problem (2.1), hence the original potential 
U alone, and even its sign, do not have a transparent 
meaning in the problem (2.12). 



C. Scattering amplitude and cross section 

Below we develop the elastic scattering theory for 
the 2D Dirac fermions in the presence of the axially- 
symmetric potential U(r). Our goal is to express the 
scattering amplitude and the cross section in terms of 
the scattering phase shifts for the spinor spherical waves 
of the form (1.20). 

First we recall that in the nonrelativistic case, with the 
incident flux along the x-direction, the 2D wave function 
has the asymptotic form [our notation follows Ref. 39] 



* ~ e lpx + 



f(0) 



x e • 



(2.16) 



B. Short-distance behavior: 
Limitations on the Dirac description 

Consider the potential U(r) that at r — > is more 
singular than 1/r. In this case, for small r, Eqs. (2.4a) 
and (2.4b) take the form 



FL - UG = and G' + UF = . 



whose solutions are 



(2.14) 



F = C sin (fUdr + S), G = C cos (f'Udr + S) 

(2-15) 

with constant C and 5. These functions strongly oscil- 
late and have no limit for r — > 0. In the nonrelativistic 
case this situation would be equivalent to falling into the 
source of the potential: namely, such a potential allows 
for infinitely deep-lying bound states. 39 Physically, a suf- 
ficiently singular potential in a (massive) relativistic sys- 
tem causes the Dirac vacuum breakdown (the Schwingcr 
effect). 40 Such a singular attractive potential will be re- 
sponsible for the free emission of electron-hole pairs; if 
the potential is attractive, electrons would then bind to 
it, while holes will be pushed to infinity. 33 If the poten- 
tial is repulsive, it will push away the electrons and bind 
holes instead. 

In the Schrodinger case falling into the potential center 
first occurs for the 1 /r 2 singularity. 39 It is not surprising 
that the Dirac problem is more sensitive to singular be- 
havior at short distances, as it can be roughly thought of 
as a "square root" of the Schrodinger equation. 

As a result of these simple considerations, for both 
repulsion and attraction, the potentials that are more 
singular than 1/r for r — > cannot be correctly consid- 
ered within the low-energy effective Dirac theory (2.1). 
In this case the exact eigenstates have to be determined 
on the length scale of the graphene lattice, where the 
long- wavelength description (2.1) breaks down. Such a 
situation, where the effect of the lattice cannot be simply 
incorporated by means of the effective-mass description, 
is reminiscent of that for the deep-lying impurity levels 
in the middle of the band gap in a semiconductor, where 
the effective- mass theory is inapplicable from the outset. 



where / is the 2D scattering amplitude, and the factor 
= e~ l7r / 4 is introduced for further convenience. The 
differential and the total cross sections, that have the 
dimensionality of length, are 39 



d 4 = \m\\ 



A= <b\f(e)\U6. 



(2.17) 



[We have denoted the scattering cross section by A since 
the letter a is commonly reserved for the conductivity] 
One way to find the scattering amplitude / is to represent 
the wave function \& in the spherical wave basis, ^ = 
^l m A m R pm {r)<& m {Q), and to consider the Schrodinger 
equation 



1 



Z {rR')' + 



' 2 m 2 2MU(r) 
P ~ r 2 ~ n 2 



R pm = (2.18) 



for each of the radial components R pm . The scattering 
phase shifts S m are then defined by the asymptotic form 
of the solutions of Eq. (2.18): 

R pm (r) ~ cos (pr - ™ - ^ + 5 m ) . (2.19) 

Using the decomposition of the plane wave 

oo 

e i P x = i m J m (pr)e imd , (2.20) 

m— — OO 

together with the definition (2.16), find 

A m = i m p~h iS ™ , (2.21) 

and 

CO 

/(#) = —= = V (s m -iy m0 , S m = e 2 ^. 

v 1 771 — — OO 

(2.22) 

From Eq. (2.22), the total cross section (2.17) 
follows: 41 ' 42 



A = - sin 2 5n 



(2.23) 



m— — oo 
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and the momentum-relaxation (transport) cross section 

2 

(2.24) 



A tr = <bd6(l-cos6)\f(6)\ 2 = - ^ sm z (S m+1 -Sr, 

J P m=-oo 

With our definition of /, the 2D optical theorem is then 



A = v^Tr/px Im/(0). 



(2.25) 



Turning to the Dirac case with the Hamiltonian (2.1), 
the asymptotic form for the spinor wave function is 



tp = u e 



+ 



f(0) 



x u f 



(2.26) 



where pe = p(cos 9, sin 9) defines the direction of scat- 
tering, and u e>p is the normalized plane wave amplitude 
(1.9). Since, according to Eq. (1.5), the scattered current 



l/l 5 



,t 



l/| 2 P* 



(2.27) 



and the incident current J in = px/e, it follows that the 
scattering amplitude / is analogous to that in the nonrel- 
ativistic case, with the cross section given by Eq. (2.17). 
Similarly, one defines the scattering phase shifts Sj via 
the asymptotic form of the radial wave functions of the 
spherical spinor (1.18), 



w > \iG(r)<S> J+1/2 {9) 



(2.28) 



[here we relabeled vp m — ► tpj using (1.21)]. The wave func- 
tions F and G are determined by Eqs. (2.4), or Eqs. (2.5). 
Their asymptotic behavior should then be compared to 
that of the free spherical spinor (1.20) with R pm (r) reg- 
ular at r = 0, 




cos pr 



JIT 



+ 6, 



Sln \pr-- + 5 3 



(2.29a) 



(2.29b) 



The spinor wave function (2.26) is represented in the ba- 
sis (1.18) as V = J2j Ajipj(r), where the coefficients Aj, 
expressed in terms of the phase shifts Sj introduced in 
Eq. (2.29), are given by the "nonrelativistic" Eq. (2.21), 
Aj = V~ x l 2 p~ x l 2 exp(iJj). The scattering amplitude, 
the cross section and the optical theorem directly follow, 
cf. Eqs. (2.22), (2.23) and (2.25) correspondingly: 



f(0) = 



1 



where Sj 
A 



A tr = 



e 2 ^ ; 
4 
P 

2 
P 



E ^--l)e 

hi +3 

E Sin2 S 3 ' 

-I +3 

E sin 2 (S j+1 - Sj) . 



i(j-l/2)9 



(2.30) 

(2.31) 
(2.32) 

(2.33) 



Finally, consider the most common massless case, 
where we derive an important property 



M = : <L 



J 3 ■ 



(2.34) 



For that we apply the symmetry (2.6) to the asymptotic 
form (2.29) and use Eqs. (1.14) [noting that Sj are de- 
fined modulo 7r, since the observable quantities are the 
S'-matrix elements (2.31)]. The property (2.34) ensures 
that backscattering vanishes in the massless limit: 



M = 0: /(tt)=0. 



(2.35) 



The absence of backscattering is a result of the destruc- 
tive interference between the time-reversed scattering 
paths. This happens since the (pseudo)helicity, the eigen- 
value of per, is asymptotically conserved during scat- 
tering off a potential that does not couple the Dirac 
points. In other words, the Dirac "spin" always remains 
in the direction of the particle's momentum. 43 The time- 
reversed backscattering paths then acquire phase differ- 
ence e l7T = —1 that corresponds to the Berry's phase 
—i<fdTip^d T ip = \§drd T 9 = tt accumulated while en- 
circling the Dirac point, tp being the spinors (1.9) with 
M = 0. 43 



D. Born approximation 

Let us now consider the potential term in the Hamil- 
tonian (2.1) as a perturbation, and find the scattering 
amplitude to the lowest order in U. Following the stan- 
dard recipe, 39 write the wave function in the form 



^ = V> (0) + v> (1) , 



(2.36) 



where the unperturbed part = e lpr u e;p is a spinor 
plane wave (1.9) in the direction of the incident momen- 
tum p, and the scattered part obeys the equation 



[W -e]V (1) = -U^ a \ 
The solution of this equation 



(2.37) 



1] = - j d2r ' Ge(r-r') \-iad T , + Ma 3 + e] U (r> £;p e lpr 

(2.38) 

is found by multiplying both sides by the operator 
Ho + e, and evaluating the operator inverse G £ = 
[Hq — (e + zO signe) 2 ] in the Fourier space [signe se- 
lects the retarded or the advanced part, corresponding to 
the particle and hole sectors]: 



G e (r) 



d 2 k 



(2tt) 2 k 2 +M 2 - (e + iQ signe) 2 



47T 



x m signe x Hq (pr) 



where p = +\ / e 2 — M 2 , and 

H^ 2 \x) = J m (x)±iY m (x) 



(2.39) 



(2.40) 
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are the Hankel's functions of the first and second kind. 
The asymptotic form of the Green's function (2.39) fol- 
lows from Eqs. (1.14), 



i / 2 

G±|e|(r) ^ ±-\ x ( 

4 V trpr 



ipr—iir J A 



(2.41) 



We now substitute the asymptotic form (2.41) into 
Eq. (2.38), applying the standard approximation |r — 
r'| r — r ■ r' and p\r — r'\ + pr' w pr — qr', where 
p' = pf is the momentum scattered in the direction of 
observation and q = p' p is the momentum transfer. 
Next we integrate by parts to switch the derivative 

e- lp ' r ' \-iad T , + Ma 3 + c] e lpr ' -» [<rp' + M<r 3 + e] e~ iciT 

use [<rp' + MiT3 + e] u±| e | ;p = ±pb(6)u±\ e \ ;p i, and com- 
pare the resulting asymptotic form to Eq. (2.26) in order 
to obtain 




Here 6 = Z(p', p), q = 2psin(6>/2), C/ q = Jdr e -^ r U(r), 
and the factor (2.43) comes from the spinor structure of 
the eigenstates (1.9). We restored h and v so that p are 
wave vectors, with (hvp) 2 = e 2 — {Mv 2 ) 2 , to make it 
explicit that / has the dimension of [length] 1 / 2 . 

In the massless limit (pristine graphcnc), the factor 
b{9) reduces to the familiar expression coming from the 
Berry phase, 43 yielding 



|- x C/ q (1 + e~ ie ) . (2.44) 



For M = the backscattcring is absent in agreement 
with the general property (2.35). 

In the nonrelativistic limit e ~ Mv 2 + (hp) 2 /2M, both 
the spinor part u p —> (1 0) T and the Berry phase factor 
b(6) — ► 2Mv/hp become trivial, yielding 39 



fT n (o) = -- 



M 



x U a 



(2.45) 









fp^)M\/M + e 



FIG. 1: Analytic continuation in the energy domain 

where Hm' 2 ^ are the Hankel's functions (2.40). The func- 
tions A{e) and B{e) become uniquely defined on the 
physical sheet of the Riemann surface of the square root 
(Fig. 1) described below. On the physical sheet, the so- 
lutions of Eqs. (2.8) can be obtained from those of (2.10) 
by analytic continuation. 

Consider the complex plane of e (Fig. 1) with the 
branch cuts along the real axis connecting the points 
e = ±M with infinity. The states on the branch cuts 
correspond to the continuous spectrum, while the real 
poles in the interval — M < e < M correspond to the 
bound states. Define \/M - e > and \[M + e > 
for —M < e < M. Analytic continuation onto e > M 
and e < — M should agree with the standard causality 
arguments (particles propagate forward in time). Since 
the time evolution ~ 9(t)e~ tept of the particle states 
(Ree > M) is described by the retarded Green's func- 
tion G R (e) ~ (e — e p + i0) _1 [here 9{t) is a unit step 
function], the branch cut to the right of M is shifted 
by the infinitesimal amount — iO below the real axis, 
Ime < 0. The square root \J e — M for e > M then 
has to be continued from the upper side of the cut, 
\fe~- ~. M — > i\/M — e, as in the Schrodinger case. 39 Con- 
versely, the hole states (Ree < —M) are governed by 
the advanced propagator G A {e) ~ (e — e p — iO)^ 1 , such 
that JdeG A (e)e- iet - 9{-t) e - l ^\ effectively shifting 
the other cut above the real axis, Ime > 0. The square 
root \f—(e + M) for e < — M then has to be continued 
from the lower side of the cut, y/—(e + M) — > i\ft + M. 
Summarizing, 



e > : 
e < : 



VMT~e 



-We - M ; (2.48a) 
-iy/-(e + M) . (2.48b) 



This determines the sign of the continuation 



E. Analytical properties 

Here we derive a few properties of the scattering solu- 
tions by considering them as functions of energy e in the 
complex plane. 39 Consider the r — > oo asymptotic form 
of the solution to the radial equations of Sec. II A, 

F ~ ,4(e)i?+ {pr) + B(e)R~ {pr) , (2.46) 
p(e) - y/(e-M)(e + M). Here [cf. Eqs. (1.14)] 

R±{pr) = ^ ±inj/2H ^' 2) iPr) * ^ ±wr , (2-47) 



A = yjM 2 -e 2 -> -ip , p = Ve 2 - M 2 . (2.48c) 

Note that the ± signs in front of the square root in 
Eq. (2.7) [that agree with those in the free spinor (1.20)] 
appear naturally as a result of the procedure (2.48). 

The radial functions of the bound states decay for 
r — > oo. This means that the discrete spectrum cor- 
responds to zeros of the function B(e) [cf. Eq. (2.47)]. 
The functions A(e) and B{e) are connected to the par- 
tial scattering amplitudes (2.31). Indeed, comparing the 
asymptotic form (2.29a) with (2.46), obtain 



A{e)/B(e) = e 2 «^)-^ 



(2.49) 
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Thus the amplitude Sj has a pole for any bound state e b . 
Following Refs. 39 and 33, we now express the residue of 
Sj in this pole via the value A{e\,). 

Consider the radial equations (2.5). Differentiating 
them with respect to e, obtain 

{d f F^F)' r - ] -d t F^F + (e + M - U)d e Gy/r = -G*V^, 
r 

(d f G^)' r + J -d £ GV^ - (c - M - U)d e Fy/¥ = Fy/r. 
r 

Multiply the first one by — 2y/rG, the second one by 
2yfFF, then multiply Eq. (2.5a) by yf?G and Eq. (2.5b) 
by —y/rF, and add up all the four equations. After many 
terms cancel, what is left can be cast in the following 
form: 



d r [r {Fd e G - Gd e F)} = r(F 2 + G 2 



(2.50) 



Next, integrate (2.50) with respect to r from r = to r, 
having in mind the limit r —> oo. The right-hand side 
becomes unity due to normalization, while in the left- 
hand side we use the asymptotic relation 



(iV?);~ -(e + M)GV^ 



(2.51) 



that follows from Eq. (2.5a) if one neglects terms with 
U(r) and j/r. The relation (2.51) allows us to rewrite 
Eq. (2.50) for r —* oo in terms of the component F only, 

(Fy/r)' r (d e Fy/f) - (F^)(deF^)' r ^e + M. (2.52) 

Now consider the asymptotic form (2.46), where we set 
A(e) ~ A(e h ) and B{e) ~ (3(e - e b ), (3 = [dB/de] e = eh . 
Substituting it into Eq. (2.52), find 



= -. 



l 



M + e b 



2A(e h ) V M - e b ' 
Using (2.49), finally obtain the S- matrix residue 



e - e b V M + e h 



(2.53) 



(2.54) 



in terms of the coefficient ^4(e b ) in the asymptotic form 
(2.46) of the wave function. We will use the result (2.54) 
in Sec. IV to normalize the bound state wave functions. 



III. LOW ENERGY SCATTERING 

As an application of the developed formalism, consider 
scattering off a potential U{r) localized within the do- 
main of the size ~ I much greater than the graphcne lat- 
tice constant. For concreteness, take U(r) — Vo9(£ — r). 
We will be primarily interested in the situation when the 
range t is small compared to the wavelength, pi < 1. 

First note that at large distances, r > £, where the po- 
tential U(r) does not contribute, the radial components 



of the spinor (1.18) are linear combinations of the corre- 
sponding free solutions (1.20) 



F\ r> e = GyJ\t + M\ {J m (pr) cos 5j - Y m (pr) sinSj} , 

(3.1a) 



G\ r> e = ±Cy/\e-M\ {J m+ i(pr) cos Sj -Y m+1 (pr) sinSj} . 

(3.1b) 

where j = m + \ and C = \fwp~J\t\. 

For r < £, the regular at r = solutions take the form 



F\ r<e - Cy/\e + M\J m (pr), (3.2a) 
G\ r<l = sign?x Cy/\e- M\J m+1 (pr) . (3.2b) 



Here we introduced e = e — Vq and p = \JV- — M 2 . The 
solutions (3.2) are written for |e] > M. For \e\ < M, 
p — > iX [cf. Sec. HE, e — > e], their counterparts read 



F\ r<e = CVM + tl m (\r), 
G\ r<e = -C'VM^ll m+1 (\r), 



(3.2c) 
(3. 2d) 



where I m (x) = i 171 J m (ix) is the modified Bcssel function 
of the first kind; C and C' are constants that can be 
chosen real. 

Both spinor components must be continuous at r = I. 
This translates into the following matching condition 



F 
G 

yielding (for \e\ > M) 



(3.3) 



7 J ™(pt) =r J m (p£)-Y m (p£)t&nSj 

^ J m+l {p£) J m+1 (p£)-Y m+1 (p£)t a nSj' { ■ ' 




( = sign e x 



M 



M 



(3.5) 



For |?| < M, C 



A/4 



M- 



[cf. Sec. HE], and 



the left-hand side of Eq. (3.4) is analytically continued 

to — \J x Im{M)/I m +\{\t), in accord with what 

one gets by applying the condition (3.3) directly to 
Eqs. (3.2c) and (3.2d). Hence we will work with Eq. (3.4) 
keeping in mind this analytic continuation for |ej < M. 
As a result, from Eq. (3.4) find [cf. Eq. (2.40)] 



5,-l = 2x 



m+l 



(J m (p£)H£Up£) - (Jm+l(pe)Hin\p£) 



(3.6) 

For the massless case [(,( — > ±], Eq. (3.6) has been first 
obtained in Ref. 26 (see also Ref. 27). 

For short-ranged scatterers, p£ -C 1, the j = ±5 
channels provide the main contribution. Indeed, for 
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m ^ 0,-1 the corresponding Sj — 1 are small as powers 
of p£; this can be seen from the short-distance behavior 
(1.15). The scattering amplitude can be thus approxi- 
mated by taking into account only the j = ±| channels, 



C Ji 7EP<? C P"i(P £ ) ~YEp£ 

(3.7) 

Here wc utilized the asymptotic behavior of the Hankel's 
functions (2.40) 



7T 



(3.8) 



In the massless limit, relevant for scattering off short- 
ranged impurities in pristine graphene, the H^\pt) con- 
tribution dominates, and the amplitude (3.7) is small as 



-1/2 



X pf(l 



,-i8\ 26 



resulting in negligible scat- 



tering away from resonance [Jo(p£) ^ 0]. The angular 
distribution of the scattered particles has a distinctive 
cos 2 (#/2) dependence that comes from the spinor struc- 
ture of the eigenstates (1.9), and agrees with the general 
property (2.35). 

In the opposite, nonrelativistic limit, C ^> (, the 
first term (m = channel) determines the amplitude 
(3.7), where now only the logarithmic term coming from 



Hq (pi) can be kept in the denominator 



39 



IV. COULOMB SCATTERING 
A. Short distance behavior. Critical field strength 

The analysis of Sec. II B shows that the low energy 
Dirac theory is inconsistent with singular potentials U ~ 
r~ s , s > 1. Hence, it is clear that the Coulomb potential 



U(r) = —hv x — 
r 



(4.1) 



is a borderline case which should be studied with care: 
Any slightly more singular potential at r — > would 
cause the Dirac vacuum breakdown. Below we consider 
the potential (4.1) where the strength a can be both pos- 
itive (attraction) and negative (repulsion). 

Suppose for now that the effective impurity strength 
a is sufficiently small, and consider Eqs. (2.5). Taking 
the short-distance behavior of the radial wave function 
as F^/r <~ r 7 and Gy/r <~ r 7 , and neglecting the non- 
singular terms as r — > 0, obtain 
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VJ 2 - a 2 



J - 1 - 2 ' 2 ' "•" ' 



(4.2) 



The theory (2.1) and (4.1) is then consistent when |a| < 
|j| for any possible j, i.e. under the condition 



a < 



(4.3) 



For imaginary 7, i.e. when \a\ exceeds \j\, the eigenstates 
F and G oscillate and have no well-defined limit as r — > 0, 
which corresponds to the Dirac vacuum breakdown in 
the same sense as in the discussion of Sec. II B. Such an 
upper bound on the potential strength is similar to that 
in the nonrelativistic collapse in the 1/r 2 potential. 39 

The condition (4.3) appears to be even more restrictive 
than that in 3D, where the Dirac theory with a point-like 
Coulomb potential source is consistent for |a| < l. 33 The 
problem of what happens when a > 1 in 3D has been 
the subject of intense theoretical investigation. 35-38 Clas- 
sically, this instability corresponds to falling of a K-shell 
electron into the potential center. On a quantum level, 
the Dirac vacuum breaks down by a sufficiently strong 
Coulomb center with Z above a certain value Z c , by cre- 
ating electron-positron pairs; an electron then binds to 
the nucleus while a positron flies off to infinity. The ma- 
jor difficulty is that the Z > Z c problem requires ultra- 
violet regularization, such as introducing the finite size 
of the nucleus. 35 However, due to very small value of the 
fine structure constant e 2 /he = 1/137, the consequences 
of this restriction never materialized in QED for the K- 
shell electrons in heavy atoms, as Z x e 2 /he < 1 for all 
the known elements in the periodic table, Z < 110. 

In a physically relevant case when the field (4.1) is due 
to a Coulomb impurity in the vicinity of the graphene 
sheet, the bare potential strength 



a = 



hv 



2e 2 

£ + 1 



(4.4) 



Here Z is the impurity valence, e is the unit charge, and 
e is the dielectric constant of a substrate. The vacuum 
value ao|z=i,e=i ~ 2.2 for v = 1 x 10 6 m/s, while for the 
Si02 substrate, ao|z=i,e=3.9 ~ 0.9. 

Electron-electron interactions result in screening which 
generally changes the shape of the potential. This 
is what usually happens in a semiconductor with a 
parabolic band, where the Coulomb potential is cut off 
on the screening length scale. In graphene, due to the 
scmimetallic electron dispersion, the screening is unusual. 
In particular, for massless Dirac fermions at half-filling, 
the linear (RPA) screening is scale invariant: 15 ' 44 it pre- 
serves the shape of the potential, and simply reduces the 
impurity strength, 



ao 



a = ao/eRPA 



by the factor 



£rpa = 1 + TZ- x 
Ahv 



q 



1 + - x — . 

2 hv 



(4.5) 



(4.6) 



Taking literally, the linear screening yields the reduction 
by the factor £rpa|s=3.9 ~ 2.4 for an impurity strength 
in the presence of the Si02 substrate, a|z=i,£=3.9 ~ 0.4. 

As one can readily see, due to the threshold (4.3), and 
a sufficiently small Fermi velocity v w c/300, the situa- 
tion in graphene is more complex than that in QED. For 
sufficiently large values of a, especially for multivalent 
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impurities, the nonlinear screening should be applied in- 
stead of the linear (RPA) response, as the latter formally 
applies only for e 2 /Tiv -C 1. In particular, a practically 
important question is whether the threshold (4.3) can 
at all be determined within the linear screening frame- 
work (4.6), i.e. whether it applies to the screened value 
(4.5). Strictly speaking, near the threshold, where the 
bare ao ~ 1, the lattice-scale physics starts playing a role, 
while the RPA dielectric constant applies in the limit of 
large distances and weak perturbations. The definition 
of the threshold as |a |/£RPA < 1/2 may probably be 
used only as an upper estimate of the threshold value. 

An even more interesting problem is the screening in 
the massive case. The linear screening, formally valid 
for e 2 /hv <C 1, is cut off beyond the (reduced) Compton 
wavelength Ac = K/Mv, hence the shape of the poten- 
tial becomes more complex. For sufficiently weak inter- 
actions, e 2 /fiv <C 1, one may argue that the screening 
can be neglected, a ~ a a , at low energies (e.g. for 
describing the bound states), since the corresponding 
Bohr radius a B = X c /{e 2 /Uv) S> A c . Taking into ac- 
count corrections in e 2 /Tiv would then amount to the 
"fine structure" of the "atomic levels" associated with 
the impurity. On the other hand, for e 2 /hv ~ 1, the 
Bohr radius and the Compton wavelength coincide; such 
a strongly-interacting "relativistic" Dirac atom will have 
deep-lying bound states. For sufficiently strong potential 
these states will reach the hole continuum (critical im- 
purity), resulting in the vacuum breakdown. In general, 
this strong-coupling problem, that requires investigation 
of the supercritical region, involves many body treatment 
that is beyond the scope of this work. In what follows we 
will assume that the condition (4.3) holds for the effec- 
tive value of impurity strength a, and consider only the 
subcritical regime. 



B. Discrete spectrum, |e| < M 

Similar to the 3D case, 33 we look for the solutions of 
Eqs. (2.5) in the form [K = v = 1] 

F = y/M+le~ p / 2 p-'- 1 / 2 F(p), (4.7a) 
G = VW~ee- p/2 p^- 1/2 G(p). (4.7b) 



Here p = 2Xr, and A = V M 2 — e 2 . After substituting 
the functions (4.7) into (2.4), obtain the equations for F 
and G: 

pF' p + h-j)F-^F-G) + ^-G = 0,(4.8a) 

pG' p + { 1 + j)G+ P -{F-G)-^-F = 0.(4.8b) 
Representing 

F = Qi + Q 2 and G = Qi-Q 2 , (4.9) 

find 

pQ'i + (7 - y ) Qi ~ (j + ^) Q2 = 0, (4.10a) 
P Q' 2 + (7 - P + y ) Q2 - [j - ^) Qi = 0, (4.10b) 



from which the equations for Qi and Q 2 are 

pQ'l + (1 + 2 7 - p)Q[ - (7 - y) Qi = , (4.11a) 



pQ'i + (1 + 2 7 - p)Q' 2 - (1 + 7 - — J Q 2 = . (4.11b) 

To derive (4.11) we used the identity 

f - M 2 a 2 /\ 2 = 7 2 - a 2 e 2 /X 2 . (4.12) 
Eqs. (4.11) are of the Kummer form, 

zT" + (c - z)T' -aJ r = 0, (4.13) 
where T is the confluent hypergeometric function 



, a z a(a + 1) z 2 

Thus the solutions of Eqs. (4.11) 



(4.14) 



Qi - C 1 F(i-ae/\,l + 2 T ,p), (4.15a) 
O2 = C 2 T{l + 1 -ael\,l + 2 1 ;p). (4.15b) 



Using J 7 (a, c; 0) = 1 and Eqs. (4.10) we find the ratio 

7 — ae/X 



= 9i _ 

° 12 ~ Ci j + Ma/X ' 
and the wave functions of the bound states 



(4.16) 



F = \fMV~e e~ p / V~ 1/2 Ci {F(j - ae/X, 1 + 2 7 ; p) 

+c 12 T{l + 1 -at/X,l + 2 T ,p)}, (4.17a) 

G = VM - e e- p / V~ 1/2 Ci {F{ n - ae/X, 1 + 2 7 ; p) 

-c 12 T{\ + 7 - ae/X, 1 + 2 7 ; p)} , (4.17b) 

where C\ is the overall normalization factor. 

Bound states occur when the functions F reduce to 
polynomials, i.e. when 



70') - 



ae r , 



A ( £ «j) 



n = 0,1,2,... for j > 0, 
n = 1,2, 3, ... for j < 0. 



From Eq. (4.18) the bound state energies follow: 34 
M sign a 



(4.18) 



4 



1 + 



7(7) = V.f - « 2 • (4.19) 



(«+7) 2 



The bound states are doubly degenerate, e n j — e n ,-j- 

The overall normalization factor C\ can be found by 
comparing the r — > 00 asymptotic behavior of the func- 
tion (4.17a) [where the leading contribution comes only 
from the first term], with Eq. (4.34) that will be obtained 
below. The asymptotic behavior of (4.17a), 



f * ( _ rCi rq + 27)yM+^ +7 - 1/2e -A, 

V ' T(l + 2 7 + n) V ' 

is found using the formula 
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H^c;z) = f ^-^(- z )-»g(a,a-c+l,-z)+'^e z z a - c g(c-a,l-a,z), (4.20) 

n -, ac a(a+l)c(c+l) ,,„,•> 

£?a,c;2 = 1 + - + , v 2 '- + ... 4.21 

1! x z 2! x z 2 

[Ref. 39, Eq. (d.14)]. As a result, the wave functions for the bound states [the upper sign corresponds to F and the 
lower one to G; e = e n j] 



(4.22) 

The functions (4.22) are normalized to f^rdr (F 2 + G 2 ) = 1. 

The size l n j of the bound state wave functions (4.22) is controlled by the parameter A = X(e n j), as 

u \ 1 V( n + i) 2 + a2 n — ; — — 2 aB 

lM 55 A^y = (Mv/h) x H 55 V(n + 7) + a x t ' flB = ^ (4 - 23) 

scaling with the "Bohr radius" that is a ratio of the reduced Compton wavelength fi/Mv and the effective fine 
structure constant e 2 /Tiv [note again that we assume weak coupling e 2 /hv <C 1]. 

C. Continuous spectrum, e > M 

The simplest way to obtain the continuous spectrum solutions in the problem (2.1) and (4.1) is to analytically 
continue the solutions (4.17a) and (4.17b) according to the procedure (2.48). This yields p — > —2ipr, and the ratio 
(4.16) becomes 

-lit- J - ia e ae 
C12 -> e 2 ^ = . ' , a e ee — . 4.24 

,7 + iMa.jp p 

The phase £j is real due to the identity (4.12). 

Consider now the |e| < M solutions (4.17). The analytic continuation |e| < M — > |e| > M to the continuous 
spectrum, using Eqs. (2.48) and (4.24), yields 

F = Vle + MI e^V 7 " 1 / 2 ^ {e^T(j - ia e , 1 + 2 7 ; -2ipr) + e~^(l + 7 - ia e , 1 + 2 7 ; -2ipr)} , (4.25a) 
G = - M \ e ipr r^ 1/2 C[ {e*?{i - ia e , 1 + 2 7 ; -2*pr) - e^^l + 7 - ia £ , 1 + 2 7 ; -2ipr)} , (4.25b) 

where C[ is some new overall normalization factor that has to be found by matching the asymptotic behavior of the 
solutions (4.25) with Eq. (2.29). For that, we first consider the asymptotic behavior of the second terms in Eqs. (4.25). 
Using the identity [see e.g. Ref. 39, Eq. (d.10)] 

T{a,c;z) = e z T(c- a,c,-z) , (4.26) 

we transform 

T{1 + 7 - ia e , 1 + 2 7 ; -2ipr) = e- 2ipr [T(-y - ia e , 1 + 2 7 ; -2ipr)]* ■ (4.27) 
As a result, obtain the normalized eigenstates for the continuous spectrum 



F = ^^ l %^ l|r( ^ ( + 1 ^ )l e-/ 2 (2 P r) 7 Re {e-+^ ( 7 - ia e , 1 + 2 7 ; -Xpr)} , (4.28a) 
G = ±-| > 2™ 1 IVj)^ e "" e/2 ( 2 ^) 7 Im {e ipr+ ^(7-^,l + 2 7 ;-2i P r)}. (4.28b) 



Note that for e < —M, the analytic continuation (2.48b) produces an extra minus sign for G (here ± = signe), as 
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expected from the asymptotic behavior (2.29). One can 
prove that the solutions (4.28) are correctly normalized 
by using the formula (4.20). The asymptotic pr 
behavior of the normalized solutions (4.28) 



oo 



F ~ ^L^J^jlp cos (pr-j7r/2 + a £ In2pr + 



G ~ ±^^^sm(pr-jw/2 + a e ln2pr + S j ) 

(4.29) 

deviates from that of Eq. (2.29) by the familiar loga- 
rithmically divergent Coulomb phase ln(2pr), ubiquitous 
in both the nonrclativistic 39 ' 41,42 and the relativistic 33 
cases. The scattering phases are then 



Sj =Zi+2 {] - 7) - arg T(l + 7 + iae/p) , (4.30) 
with the corresponding 5-matrix elements (2.22) 

s , = e 2tS :j = J + iMa/p r(l + 7 - iae/p) 
7 — iae/p r(l + 7 + iae/p) 

(4.31) 

As expected, the poles of Sj determined by the Gamma- 
function in the numerator of Eq. (4.31) for 1+7 — iae/p = 
1 — n, n = 1, 2, as well as by 7 — iae/p = for j > 0, 
that occur for the imaginary p = iX, give the correspond- 
ing bound states (4.19). The residues at these poles are 



in+1 



X 3 (j + Ma/\)e i <i-^ 
aM 2 n\ T(l + 27 + n)(e - e n>j ) 



(4.32) 



We now derive the asymptotic form of the discrete 
spectrum wave functions based on the relation (2.54) 
between the S-matrix residue and the coefficient A in 
the asymptotic form of the wave function (2.46). In 
the case of the Coulomb scattering, the coefficient A 
will itself depend on r due to the logarithmically di- 
vergent Coulomb phase. In particular, the left-hand 
side of Eq. (2.54) should be corrected by the factor 
e 2 M£ i„2 P r _^ (-)« e ^7(2Ar) 2 ("+T). As a result, near the 
pole 



Sje 2ia ' ln2pr — ► — e" 1 - 7 '- 



A 3 (j + Ma/A)(2Ar) 2 ("+T) 



aM 2 n\T{\ + 2 1 + n){e-e n . J ) ' 

(4.33) 

which in turn equals the right-hand side of Eq. (2.54). 
Thus we obtain the asymptotic form 



(4.34) 



A(r) 



h\ (Mv 2 + e){j + aMv/hX) 
MwV 2hvan\T{\ + 27 + n) 



x (2Ar) 



n+7 



Here we restored h and v, with HX/v — V 1 M 2 — e 2 , so 
that the dimension of F is explicitly 1/ [length]. Note 
also that \j\ < aMv/hX(e n _j), so that A(r) is always 
real. 



1. Nonrelativistic limit (parabolic band) 

The nonrelativistic limit occurs when the "nonrela- 
tivistic velocity" w nr of the particle is much smaller than 
the graphene Fermi velocity ("speed of light" v) [which 
we write here explicitly] , 



hp 



e-Mv 2 ~ M^«M V 2 . 



M 2M 

(4.35) 

In the limit v — > 00, the fine structure constant a — > 0, 
whereas the "nonrelativistic fine structure constant" 



a Dr = 



(4.36) 



remains finite, and determines the Coulomb interaction 
strength. In this case 



fipv 



Mv 
Tip 



7 -> = RI + 5 sign j , 
(4.37) 

and using e %lz ^~^ = sign j and T(l + z) = zT(z), obtain 

(4.38) 



q nr _ e 2i6™ _ + 2 ~ 



r(|m| + 5 + ZQfnr) 

It is instructive to show by a direct calculation, given 
below, that the phases (4.38) agree with the asymptotic 
behavior of the radial Schrodinger wave function R m (r) 
of the corresponding nonrelativistic Coulomb problem. 
The radial equation in the presence of the Coulomb po- 
tential U(r) = —Ze 2 /r reads 



1 d f r±R 



m" „ 2MZel „ 2Me n 



-R 



-R 



■R. (4.39) 



r dr \ dr J r' z h z r Ti z 

We begin with the bound states. After the substitution 



R = p\™\ e -i>/ 2 Q(p), p = 2A nr r, HX m = y/-2Me ni 

(4.40) 

the problem is reduced to the Kummer equation 

P Q" + (2\m\ + 1 - p)Q' - (\m\ + \ - a)Q = , (4.41) 
where 

Z MZe 2 



a = a(e m ) 



(4.42) 



a_BA nr h A nr 
Its solutions 

R = C nr x (2A nr r)l m l e - A '» r F(-n, 2\m\ + 1; 2A nr r), 
-n=\m\ + \- 5(e nr ) , (4.43) 
r _ ( xn V2A nr V(2|m|+n)! 

(2|m|)! v /n! (n + |m| + |) 

become normalizable, J^rdr R 2 (r) = 1, with 
J-(-n,2|m| + l;2A nr r) = J 2 '^'?' x L^"' ,(2A nr r) 



[(n + 2|m|)!] 2 n + 2 IH v 



(4.44) 
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given by the associated Laguerre polynomials, when n = 
0, 1, 2, yielding the nonrelativistic spectrum 45 



Z 2 Mei 



1 



,2 ' 



n = 0,1,2,... 



2^ (n+\ m \+iy 

(4.45) 

in agreement with the corresponding limit of the bound 
states (4.19). [Taking into account the discarded terms 
~ e 2 /Tiv would correspond to the "fine structure" of the 
energy levels (4.45).] 

The normalization coefficient C nr for the bound state 
solutions (4.43) is found, as above, via comparing their 
r — > oo asymptotic behavior with the nonrelativistic limit 
of Eq. (4.34), 



R 



V / 2A nr (2A„ r r)"+l m l 



-Anr^ 



r— >oo — 



. (4.46) 



^Jnl (n + 2\m\)\(n+ \m\ + \) 

While obtaining Eq. (4.46) we identified [cf. Eq. (4.18)] 

aMv ae„ 



h\ ' UvX 



a = n + \j\ = n - 



(4.47) 



The analytic continuation of the solutions (4.43) via 
A nr — > —ip, a — > janr, and subsequent asymptotic ex- 
pansion using both terms in the formula (4.20), yields 



1 ( i r> rnr TO7F 77 \ 

R m ~ — cos (pr + a nr In 2pr + d TO -J, 

(4.48) 

where the scattering phase shifts [defined mod n] 



C = - ar g r(|m| + \ + ia nl ) (4.49) 
correspond to the S'-matrix elements (4.38). 

2. Ultrarelativistic limit (graphene) 

In the ultrarelativistic limit |e| ^> Mv 2 relevant for 
pristine graphene monolayer [M = 0], 



a £ — ae/pv — > a sign e . 



(4.50) 



In this case, similar to the 3D Dirac fermions, 33 the S'- 
matrix elements (4.31) become independent of the abso- 
lute value of energy (depending only on its sign): 

e 2is j= 3 r(l+ 7 -ia e ) eM ._ 7) 

7- ia e r(l +7 + ia e ) J J 

(4.51) 

Note that the general property (2.34) for masslcss 
fermions holds. 



D. Scattering cross section 

We start from the nonrelativistic limit, where one can 
directly sum the series (2.22) with S m from Eq. (4.38) 



to obtain the nonrelativistic scattering amplitude in the 
closed form 41 ' 42 



f(0) = 



r(i - ia Dr ) 



,ia nr lnsin 2 (e/2) 



T(ia nr ) 



^2p sin 2 (6/2) 

(4-52) 

For completeness, the details are given in the Ap- 
pendix A. Using the property T(z)T(l — z) = n/ sin(7rz), 
the 2D Rutherford cross section follows 41,42 



a nr tanh ira v 



d6 



2pS\TL 



2 9 



(4.53) 



Here 9 is the scattering angle, and the momentum trans- 
fer q = 2psin|. The cross section (4.53) is written in 
terms of the nonrelativistic fine structure constant (4.36). 

In the opposite, massless limit, the phases (4.51) be- 
come independent on the absolute magnitude |e| of en- 
ergy. Thus the differential scattering cross section scales 
with the particle wavelength, 



dA 
d6 



-(&) 



(4.54) 



|e|»M 



where r(9) is an |e|-independent function of the scattering 
angle. As the symmetry (2.34) is fulfilled, the backscat- 
tering is absent: f(ir) = and t(tt) = 0, cf. Eq. (2.35). 

For the general case, the differential cross section is 
obtained by summing the series (2.30) with Sj from 
Eq. (4.31). This problem is notoriously cumbersome, as 
it has long been known from the three dimensions. 33 Un- 
fortunately, for the full relativistic problem (Mv 2 ^ oo), 
neither the differential, nor the transport cross section 
can be obtained in the closed form. Moreover, the series 

(2.30) for the total cross section with the phase shifts 

(4.31) does not converge. To obtain the converging ex- 
pression for the scattering amplitude, one needs to per- 
form an appropriate resummation of this series. 33,46 In 
Appendix B we show how to represent the Rutherford 
scattering amplitude via the convergent double integral, 
where the energy and angular dependences are separated. 
In the following section we numerically sum the series 
(2.33) for the transport cross section. 



V. EXACT COULOMB TRANSPORT CROSS 
SECTION IN GRAPHENE 

For a half-filled 7r-electron band, the RPA screening 
(4.5) is scale-invariant, preserving the functional form of 
the potential. This is why the exact solution for the 
Coulomb potential (4.1) can be practically important, as 
it survives the interaction effects at least on the linear 
screening level. 

The knowledge of the scattering phases (4.51) allows 
us to obtain the exact transport cross section (2.33) for 
scattering off the Coulomb impurity located in the im- 
mediate vicinity of the graphene sheet. Since the phase 
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shifts are energy independent for M = 0, the transport 
cross section is proportional to the energy-dependent car- 
rier wavelength A e , 

A tr = C(a e ) x A e , A e = 2nhv/\e\ . (5.1) 

The dimensionless function C(a t ) [a e = a sign e], which 
is the transport cross section in the units of the carrier 
wavelength, is plotted in Fig. 2. 

The transport cross section (5.1) has a few distinct 
features: (i) it is not symmetric with respect to the 
sign of the potential as seen by the carrier: A positively 
charged impurity (a > 0) scatters conduction electrons 
(e > 0) more effectively than it scatters holes (e < 0). 
The scattering asymmetry with respect to the sign of 
the potential arises naturally [e.g. in the next-to-leading 
Born approximation for low-enegy particles, Problem 6 
in Sec. 132 of Ref. 39]. Physically, one may expect the 
particle to spend more time around an attractive poten- 
tial center and thereby be more significantly deflected 
(although for an ultrarelativistic particle this intuition 
may fail). However, for the practically important 2D 
and 3D Coulomb scattering in a parabolic band, the cor- 
responding exact solutions are somewhat exceptional in 
a sense that they lack such an asymmetry. Remarkably, 
for the "relativistic" carrier dispersion, characteristic of 
graphene, this generally expected asymmetry is recov- 
ered, (ii) The cross section is apparently non-monotonic 
for the attractive Coulomb scatterers. (iii) The unitary 
limit, Sijt—i/2 = 7r/2 for the j = ±| partial wave, is 
reached for the mutual attraction when 

a* =a X signe » 0.494, (5.2) 

just below criticality. 

We were not able to link the above unitarity to any 
resonance or other special behavior at the point (5.2), 
which may as well be accidental. (One may argue that 
after subtracting the logarithmically divergent Coulomb 
phase, the phase shifts alone have lost their meaning, 
whereas the differences between them correspond to ob- 
servable quantities.) Indeed, the partial Coulomb scat- 
tering phases Sj even in the nonrelativistic Rutherford 
problem 39 generally pass the value it/ 2 at non-special 
values of parameters. At that point, the particular angu- 
lar momentum channel reaches unitarity (maximum pos- 
sible scattering). However, in previously studied cases 
this unitarity in one channel did not cause a local maxi- 
mum for the sum (2.24) over all channels. The opposite 
situation apparently happens in the 2D Dirac case: the 
relatively strong dependence of the scattering cross sec- 
tion on the lowest-j phase shift causes the local maximum 
shown in Fig. 2 inset. 

The conductivity of graphene monolayer in the pres- 
ence of charged impurities with the transport cross sec- 
tion (5.1) is obtained in Ref. 48. The attraction-repulsion 
asymmetry of the cross section can in principle allow one 
to determine the numbers of positively and negatively 
charged impurities independently. 



0.8 r 




-1.5 -0.4 -0.3 -0.2 -0.1 0.1 0.2 0.3 0.4 0.5 
repulsion Impurity charge a £ = a x sign e attraction 



FIG. 2: (Color online) Transport cross section as a function 
of the impurity charge a (solid line). Dashed line is the Born 
approximation; thin red line is the lowest scattering phase 
shift 8 ±1/2 /ir for |j| = 1/2. Note that S ±1/2 = tt/2 (the 
unitary limit) for a signe = a* w 0.494; the transport cross 
section decreases for a* < a e < 1/2, as shown in the inset. 
Attraction here means mutual attraction between the charged 
carrier and the impurity, i.e. sign a = signe, while repulsion 
stands for sign a = — sign e. 

Born approximation 

We would now like to compare the exact re- 
sult obtained above with the previously used Born 
approximation. 15-17 The Born scattering amplitude 
(2.44) is straightforwardly found using {7 q = —2-Khva/q: 

Thus the differential cross section 

= x cot z - , (5.4) 

dd 2p 2 ' y ' 

and the transport cross section 

2 2 

A Bor„ = = C Born (a) x ^ > ^Born = 1^2 (5 5) 

Note that the differential cross section is singular for 9 = 
as is expected from the long-range character of the 
Coulomb field. 

The above cross sections can also be obtained from the 
Golden Rule. For completeness, we present such a calcu- 
lation for the momentum relaxation time in the presence 
of Ui Coulomb impurities per unit area: 

r^HT) = ^J0^\M pp ^S(e p ~e p ,)(l~cose) 
= niir 2 (Zel) 2 /e — > niir 2 (hva) 2 /e , (5-6) 
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where the intra-band matrix element of the interaction 
taken between the spinor plane wave states (1.9) 

M pp > = J dr^. p ,*7(r)V e;p = \ (l + e 40 ) t/ q . (5.7) 

Here the Fourier transform of the effective potential (4.1) 
is Uq — —2TrZe 2 /q, and in Eq. (5.6) we accounted for 
RPA screening via the procedure (4.5). The transport 
cross section (5.5) follows from niWTt r A tr = 1. We remark 
in passing on the following curious observation specific for 
the Born approximation: The ultrarelativistic Coulomb 
transport rate (5.6) written in terms of the quasiparticle 
energy e formally coincides with the corresponding non- 
relativistic value in two dimensions, where e = p 2 /2m*. 

We also note that, as it is generally expected for the po- 
tential scattering, 39 the Born approximation (5.5) over- 
estimates the exact result for the repulsion and underes- 
timates it for the attraction. Fig. 2 shows that, numeri- 
cally, Born approximation works well for a < 0.1, while 
for the experimentally relevant values a ~ 0.5 it fails 
by about a factor of two. On a qualitative level, since 
the Born scattering assumes small phase shifts, it fails 
to recognize the strong repulsion/ attraction asymmetry, 
the unitary scattering occurring at the value (5.2), and 
the associated non-monotonic dependence of the cross 
section for the mutually attracting carrier and impurity. 

SUMMARY 

In this work we have outlined in detail the elastic scat- 
tering theory for the (2+l)-dimensional massive Dirac 
fcrmions in an axially-symmetric potential. The formal- 
ism is relevant for the transport in pristine graphene 
monolayers (massless limit), and for graphene layers with 
the broken symmetry between the sublattices (result- 
ing in the finite Dirac mass gap), in the presence of a 
smooth potential disorder. We showed that the Dirac 
theory becomes sensitive to the lattice scale for the po- 



tentials that are more singular than 1/r as r — > 0. We 
also considered scattering off a localized potential whose 
size is smaller than the Dirac fermion wavelength (but 
larger than the graphene lattice scale) . For the Coulomb 
scattering, U = —~hva/r, the exact solution is found be- 
low the threshold \a\ < 1/2; from the scattering phase 
shifts we obtain the exact transport cross section. The 
transport cross section for the massless case (pristine 
graphene) is shown to exhibit a pronounced asymmetry 
with respect to attraction versus repulsion between the 
charge carrier and the Coulomb impurity. 

Note added: (i) Upon completion of this work we 
learned about the preprint 49 where, for the Coulomb po- 
tential, the opposite, supercritical situation \a\ > 1/2 
is discussed for the massless limit. The perturbative 
rcnormalization group treatment of Ref. 49, based on the 
single-particle Friedel sum rule, applies in the weakly- 
interacting limit e 2 /Tiv -C 1. The strong coupling limit 
of the problem, valid for e 2 /Tiv ~ 1 and large impurity 
charge Z» 1, was subsequently considered in Ref. 50, 
yielding a qualitative change in the screened potential 
profile. (ii) Also, recently, angular resolved photoe- 
mission spectroscopy measurement 51 became available, 
according to which the Dirac gap A <~ 0.1 eV opens 
up for graphene on SiC substrate. In the presence of 
Coulomb impurities this would lead to the subgap states 
(4.19) with the wave function spread over ~ 10 nm [cf. 
Eq. (4.23)]. Localized states on this length scale can be 
detected using scanning techniques. 
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APPENDIX A: NONRELATIVISTIC COULOMB SCATTERING AMPLITUDE 



Here we show how to sum the series (2.22) with S m from Eq. (4.38). For that, first use the integral representation 

T(m+\-~ia m ) B(m + \ - ia nr , 2ia nr ) 
l 

2 

Noting that S™ r m = S™, we then sum the two similar looking geometric series. Defining z = e 40 , find 



r ; TO +^-^) = B(m+j-za 2 Wnr ) = 1_ /a^i/^^^-i, (Al) 
T{m+\ + ia m ) T(2ia nr ) T(2ia nr ) J V ' K ' 



oo oo p/ 1 ■ \ ) 1 

V S m z m + V S m z- m - )* , " ^(1, \ - ia nl , \ + za nr ; z) + )f , " -^(1, § - ia nr , § + za nr ; z^) . (A2) 
n i r(o + «a nr ) T(f + ia nT ) z 

Here we used the integral representation [Ref. 47, Eq. (9.111)] 

H», b, c; z) = B{b l_ b) J/tt^i 1 ty- b ~\l zty a (A3) 
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for the hypergeometric function 



!F{a, b, c; z) = 1 



ab z a(a + 1)6(6+ 1) z 2 



1! 



c(c+l) 



2! 



(A4) 



Now use the analytic continuation of the hypergeometric series for \z\ > 1, Eq. (9.132.2) of Rcf. 47 [Eq. (e.6) of 
Rcf. 39], 

F(a,b,c;z) = T ^ T J^' a \ (-z)- a T(a,a+l-c,a+l-b-z- 1 ) + T ^^ (A5) 



r(6)r(c-o) 



T(a)T{c-b) 



to transform the second term of Eq. (A2) to become the function of z rather than z 1 . Application of the first term 
of the formula (A5) cancels the first term of the sum (A2), while the second term of Eq. (A5) yields 



E 5 - 



z =- 



r( 2 +ia nr )r( 2 ia nr ) t _._i,_i/2+» anr ^iz _ 



r(2«a nr ) 

Finally, using Eq. (9.131) of Rcf. 47 [Eq. (c.4) of Rcf. 39], 

F{a, 6, c; z) = (1 - zf- a - b T{c - a,c - b,c; z) , 
as well as the doubling formula [Ref. 47, Eq. (8.335.1)] 



^"(2 ^nrj 1 2za nr , ^ ^nr! ^-) 



2x-l 



2 2. 



-T{x)Y{x+\), 



T{2x) = 



and —(1 — z) 2 / z = 4sin 2 (0/2), one obtains the scattering amplitude (4.52). 



APPENDIX B: RELATIVISTIC COULOMB SCATTERING AMPLITUDE 



(A6) 



(A7) 



(A8) 



Consider the following transformations of the series for the scattering amplitude (2.30), z = e l6 : 



E 



m— — 00 



-1/2 



E 



7 = ±i ±^ 

.7 - 1 - 2 ' 2 ' " ' 

z -l/2 



T(l + 2za e ) y 

-,-1/2 



7) 



r(l + 2ia e ) 



r(l+ 7 + ia e ) 

/ dfct- <0 "(l-t) 2ia " /^7rtan(7r K )( K + iM)^ K r- 1 e l7r( ^ 
J Q at 



c , 2iri 7(k) 



(Bla) 

(Bib) 
(Blc) 



r 



Here M = aM/p, and in the first line we canceled the 
denominator -f—ia £ by utilizing the property T(x + 1) = 
aT(a;). 

In the next line, Eq. (Bib), we utilize the integral rep- 
resentation for the Euler's Beta function 

B(P, 9) = = jf - ^ 1 • (B2) 

The summation is reduced to the contour integration by 
virtue of the fact that n tan ttk has residues at the points 
re = ±1/2, ±3/2, ... with a value —1. The re-intcgration 
is along the contour C shown in Fig. 3. 



Finally, in the line (Blc) we have performed the inte- 
gration by parts in the t-variable, f _1 = d(f)/-f. This 
way all the remaining singular behavior in the re-plane 
[besides the residues at re — ±1/2, ±3/2, ...] is reduced 
to that of 7(k) in the denominator. The phase of the 
square root in j(k) = V 're 2 — a 2 is defined in a standard 
way, Fig. 4. The exponential function e 47r ' K ~ 7 ^ does not 
diverge for large imaginary re, which allows us to deform 
the integration contour C in Eq. (Bl) to C around the cut 
between re = ±a shown in Fig. 4. Clearly, this procedure 
is valid only if the condition (4.3) is satisfied. 
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